SPDE Toy Example
Toy dataset from Blangiardo and Cameletti (2015).
# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy,
pch = 19, asp = 1, main = 'Toy Data')

# Create a mesh for the SPDE method and then plot it.
toymesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]),
max.edge = c(0.1, 0.2))
plot(toymesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5),
pch = 19)

# Initialize SPDE projector matrix.
toyA <- inla.spde.make.A(toymesh, as.matrix(SPDEtoy[,c('s1', 's2')]))
# Initialize exponential covariance structure for SPDE.
toyspde <- inla.spde2.matern(toymesh, alpha = 2)
toy_formula <- y ~ -1 + intercept + f(spatial_field, model = toyspde)
toy_data <- list(
y = SPDEtoy$y,
intercept = rep(1, nrow(SPDEtoy)),
spatial_field = seq_len(nrow(SPDEtoy))
)
# Fit the model with INLA.
# Can't find the intercept??
toyfit <- inla(
toy_formula,
# y ~ f(spatial_field, model = toyspde),
data = toy_data,
# control.predictor = list(A = toyA, compute = TRUE)
# control.predictor = list(compute = TRUE)
)
Bei Dataset
Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005).
# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')

bei_corners <- vertices.owin(Window(bei))
bei_domain <- cbind(bei_corners$x, bei_corners$y)
bei_full_mesh <- inla.mesh.2d(cbind(bei$x, bei$y),
cutoff = 50, max.edge = c(50, 100),
loc.domain = bei_domain)
plot(bei_full_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei, pch = '.', col = 'blue')

bei_full_spdf <- as.SpatialPoints.ppp(bei)
# CHECK PRIORS!
matern_full <- inla.spde2.pcmatern(bei_full_mesh,
prior.sigma = c(0.1, 0.99),
prior.range = c(5, 0.01))
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = matern_full) + Intercept
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf)
lambda_full <- predict(bei_full_lgcp, pixels(bei_full_mesh), ~ exp(mySmooth + Intercept))
# Plot posterior means and posterior sd.
plot(lambda_full)
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(lambda_full['sd'])
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
n_quads <- 10
botleft <- cbind(runif(n_quads, 0, 950), runif(n_quads, 0, 450))
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
cbind(
botleft[r, 1] + c(0, 0, 50, 50),
botleft[r, 2] + c(0, 50, 50, 0)
)
)})
bei_win <- do.call(
union.owin,
apply(botleft, 1, function(x){return(
owin(x[1] + c(0, 50), x[2] + c(0, 50))
)})
)
bei_hole <- bei[complement.owin(bei_win)]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)
plot(bei_hole, main = 'Observed Subregion', pch = '.', cols = 'black')

bei_hole_mesh <- inla.mesh.2d(cbind(bei_hole$x, bei_hole$y),
cutoff = 50, max.edge = c(50, 100),
loc.domain = bei_domain)
plot(bei_hole_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei_hole, pch = '.', col = 'blue')

bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
# CHECK PRIORS!
matern_hole <- inla.spde2.pcmatern(bei_hole_mesh,
prior.sigma = c(0.1, 0.99),
prior.range = c(5, 0.01))
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = matern_hole) + Intercept
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf)
lambda_hole <- predict(bei_hole_lgcp, pixels(bei_hole_mesh), ~ exp(mySmooth + Intercept))
# Plot posterior means and posterior sd.
plot(lambda_hole)
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(lambda_hole['sd'])
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)

bei_samp_mesh <- inla.mesh.2d(cbind(bei_samp$x, bei_samp$y),
cutoff = 50, max.edge = c(50, 100),
loc.domain = bei_domain)
plot(bei_samp_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
plot(Window(bei_samp), border = 'blue', add = TRUE)
points(bei_samp, pch = '.', col = 'blue')

bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
# CHECK PRIORS!
matern_samp <- inla.spde2.pcmatern(bei_samp_mesh,
prior.sigma = c(0.1, 0.99),
prior.range = c(5, 0.01))
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = matern_samp) + Intercept
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf)
lambda_samp <- predict(bei_samp_lgcp, pixels(bei_samp_mesh), ~ exp(mySmooth + Intercept))
# Plot posterior means and posterior sd.
plot(lambda_samp)
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(lambda_samp['sd'])
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

LS0tCnRpdGxlOiAiU3BhdGlhbCBQcmVkaWN0aW9uIHdpdGggSU5MQSIKYXV0aG9yOiAiS2VubmV0aCBBLiBGbGFnZyIKYmlibGlvZ3JhcGh5OiAiLi4vcmVmZXJlbmNlcy5iaWIiCm91dHB1dDoKICBodG1sX25vdGVib29rOgogICAgZmlnX2hlaWdodDogNgogICAgZmlnX3dpZHRoOiAxMAogICAgZmlnX2Nyb3A6IEZBTFNFCiAgICBoZWlnaHQ6ICI5NjBweCIKICAgIHdpZHRoOiAiNzIwcHgiCiAgICBzZWxmX2NvbnRhaW5lZDogVFJVRQotLS0KCgpgYGB7ciBzZXR1cCwgY2FjaGUgPSBGQUxTRSwgZWNobyA9IEZBTFNFLCBtZXNzYWdlID0gRkFMU0UsIHdhcm5pbmcgPSBGQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGNhY2hlID0gRkFMU0UsIGVjaG8gPSBUUlVFLCB3YXJuaW5nID0gRkFMU0UsCiAgbWVzc2FnZSA9IEZBTFNFLCBkcGkgPSAxNTAsIGZpZy5hbGlnbiA9ICdjZW50ZXInKQpgYGAKCmBgYHtyIHBhY2thZ2VzLCBjYWNoZSA9IEZBTFNFLCBlY2hvID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSwgd2FybmluZyA9IEZBTFNFfQpsaWJyYXJ5KHNwYXRzdGF0KQpsaWJyYXJ5KElOTEEpCmxpYnJhcnkoaW5sYWJydSkKbGlicmFyeShtYXB0b29scykKYGBgCgoKIyBTUERFIFRveSBFeGFtcGxlCgpUb3kgZGF0YXNldCBmcm9tIEByaW5sYS4KCmBgYHtyIHNwZGV0b3ksIGZpZy53aWR0aCA9IDYsIG91dC53aWR0aCA9ICc2MCUnfQojIFBsb3QgdGhlIGRhdGEuCnBsb3QoczIgfiBzMSwgY29sID0gcmdiKFNQREV0b3kkeSAvIG1heChTUERFdG95JHkpLCAwLCAwKSwgZGF0YSA9IFNQREV0b3ksCiAgICAgcGNoID0gMTksIGFzcCA9IDEsIG1haW4gPSAnVG95IERhdGEnKQpgYGAKCmBgYHtyIHNwZGVtZXNoLCBmaWcud2lkdGggPSA2LCBvdXQud2lkdGggPSAnNjAlJ30KIyBDcmVhdGUgYSBtZXNoIGZvciB0aGUgU1BERSBtZXRob2QgYW5kIHRoZW4gcGxvdCBpdC4KdG95bWVzaCA8LSBpbmxhLm1lc2guMmQoYXMubWF0cml4KFNQREV0b3lbLGMoJ3MxJywgJ3MyJyldKSwKICAgICAgICAgICAgICAgICAgICAgICAgbWF4LmVkZ2UgPSBjKDAuMSwgMC4yKSkKcGxvdCh0b3ltZXNoLCBhc3AgPSAxKQpwb2ludHMoU1BERXRveSRzMSwgU1BERXRveSRzMiwgY29sID0gcmdiKFNQREV0b3kkeSAvIG1heChTUERFdG95JHkpLCAwLCAwLCAwLjUpLAogICAgICAgcGNoID0gMTkpCmBgYAoKYGBge3Igc3BkZWZpdH0KIyBJbml0aWFsaXplIFNQREUgcHJvamVjdG9yIG1hdHJpeC4KdG95QSA8LSBpbmxhLnNwZGUubWFrZS5BKHRveW1lc2gsIGFzLm1hdHJpeChTUERFdG95WyxjKCdzMScsICdzMicpXSkpCgojIEluaXRpYWxpemUgZXhwb25lbnRpYWwgY292YXJpYW5jZSBzdHJ1Y3R1cmUgZm9yIFNQREUuCnRveXNwZGUgPC0gaW5sYS5zcGRlMi5tYXRlcm4odG95bWVzaCwgYWxwaGEgPSAyKQoKdG95X2Zvcm11bGEgPC0gIHkgfiAtMSArIGludGVyY2VwdCArIGYoc3BhdGlhbF9maWVsZCwgbW9kZWwgPSB0b3lzcGRlKQp0b3lfZGF0YSA8LSBsaXN0KAogICAgeSA9IFNQREV0b3kkeSwKICAgIGludGVyY2VwdCA9IHJlcCgxLCBucm93KFNQREV0b3kpKSwKICAgIHNwYXRpYWxfZmllbGQgPSBzZXFfbGVuKG5yb3coU1BERXRveSkpCiAgKQoKIyBGaXQgdGhlIG1vZGVsIHdpdGggSU5MQS4KIyBDYW4ndCBmaW5kIHRoZSBpbnRlcmNlcHQ/Pwp0b3lmaXQgPC0gaW5sYSgKICB0b3lfZm9ybXVsYSwKIyAgeSB+IGYoc3BhdGlhbF9maWVsZCwgbW9kZWwgPSB0b3lzcGRlKSwKICBkYXRhID0gdG95X2RhdGEsCiMgIGNvbnRyb2wucHJlZGljdG9yID0gbGlzdChBID0gdG95QSwgY29tcHV0ZSA9IFRSVUUpCiMgIGNvbnRyb2wucHJlZGljdG9yID0gbGlzdChjb21wdXRlID0gVFJVRSkKKQpgYGAKCgojIEJlaSBEYXRhc2V0CgpFeGFtcGxlIGZyb20gQG1vZWxsZXJ3YWFnZXBldGVyc2VuLCBfQmVpbHNjaG1pZWRpYSBwZW5kdWxhIExhdXJhY2VhZV8gbG9jYXRpb25zCmluIGEgcGxvdCBpbiBQYW5hbWEuIGBiZWlgIGRhdGFzZXQgaW4gYHNwYXRzdGF0YCBbQHNwYXRzdGF0XS4KCmBgYHtyIGJlaXB0c30KIyBQbG90IHRoZSBmdWxsIHBvaW50IHBhdHRlcm4uCnBsb3QoYmVpLCBwY2ggPSAnLicsIGNvbHMgPSAnYmxhY2snLCBtYWluID0gJ1JlYWxpemVkIFBvaW50IFBhdHRlcm4nKQpgYGAKCmBgYHtyIGJlaW1lc2h9CmJlaV9jb3JuZXJzIDwtIHZlcnRpY2VzLm93aW4oV2luZG93KGJlaSkpCmJlaV9kb21haW4gPC0gY2JpbmQoYmVpX2Nvcm5lcnMkeCwgYmVpX2Nvcm5lcnMkeSkKYmVpX2Z1bGxfbWVzaCA8LSBpbmxhLm1lc2guMmQoY2JpbmQoYmVpJHgsIGJlaSR5KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgY3V0b2ZmID0gNTAsIG1heC5lZGdlID0gYyg1MCwgMTAwKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbG9jLmRvbWFpbiA9IGJlaV9kb21haW4pCnBsb3QoYmVpX2Z1bGxfbWVzaCwgYXNwID0gMSkKcGxvdChXaW5kb3coYmVpKSwgYm9yZGVyID0gJ2JsdWUnLCBhZGQgPSBUUlVFKQpwb2ludHMoYmVpLCBwY2ggPSAnLicsIGNvbCA9ICdibHVlJykKYGBgCgpgYGB7ciBiZWlmdWxsbGdjcCwgY2FjaGUgPSBUUlVFfQpiZWlfZnVsbF9zcGRmIDwtIGFzLlNwYXRpYWxQb2ludHMucHBwKGJlaSkKIyBDSEVDSyBQUklPUlMhCm1hdGVybl9mdWxsIDwtIGlubGEuc3BkZTIucGNtYXRlcm4oYmVpX2Z1bGxfbWVzaCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwcmlvci5zaWdtYSA9IGMoMC4xLCAwLjk5KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwcmlvci5yYW5nZSA9IGMoNSwgMC4wMSkpCmNtcF9mdWxsIDwtIGNvb3JkaW5hdGVzIH4gbXlTbW9vdGgobWFwID0gY29vcmRpbmF0ZXMsIG1vZGVsID0gbWF0ZXJuX2Z1bGwpICsgSW50ZXJjZXB0CmJlaV9mdWxsX2xnY3AgPC0gbGdjcChjbXBfZnVsbCwgYmVpX2Z1bGxfc3BkZikKbGFtYmRhX2Z1bGwgPC0gcHJlZGljdChiZWlfZnVsbF9sZ2NwLCBwaXhlbHMoYmVpX2Z1bGxfbWVzaCksIH4gZXhwKG15U21vb3RoICsgSW50ZXJjZXB0KSkKCiMgUGxvdCBwb3N0ZXJpb3IgbWVhbnMgYW5kIHBvc3RlcmlvciBzZC4KcGxvdChsYW1iZGFfZnVsbCkKcGxvdChXaW5kb3coYmVpKSwgYm9yZGVyID0gJ3doaXRlJywgYWRkID0gVFJVRSkKcG9pbnRzKGJlaSwgcGNoID0gJy4nLCBjb2wgPSAnd2hpdGUnKQpwbG90KGxhbWJkYV9mdWxsWydzZCddKQpwbG90KFdpbmRvdyhiZWkpLCBib3JkZXIgPSAnd2hpdGUnLCBhZGQgPSBUUlVFKQpwb2ludHMoYmVpLCBwY2ggPSAnLicsIGNvbCA9ICd3aGl0ZScpCmBgYAoKYGBge3IgYmVpaG9sZX0KIyBUYWtlIGEgc2FtcGxlIG9mIHF1YWRyYXRzIGFuZCBwbG90IHRoZSBvYnNlcnZlZCBwb2ludCBwYXR0ZXJuLgpzZXQuc2VlZCg4NDMyMykKbl9xdWFkcyA8LSAxMApib3RsZWZ0IDwtIGNiaW5kKHJ1bmlmKG5fcXVhZHMsIDAsIDk1MCksIHJ1bmlmKG5fcXVhZHMsIDAsIDQ1MCkpCmJlaV9pbnRlcmlvciA8LSBsYXBwbHkoc2VxX2xlbihucm93KGJvdGxlZnQpKSwgZnVuY3Rpb24ocil7cmV0dXJuKAogICAgY2JpbmQoCiAgICAgIGJvdGxlZnRbciwgMV0gKyBjKDAsIDAsIDUwLCA1MCksCiAgICAgIGJvdGxlZnRbciwgMl0gKyBjKDAsIDUwLCA1MCwgMCkKICAgICkKICApfSkKYmVpX3dpbiA8LSBkby5jYWxsKAogIHVuaW9uLm93aW4sCiAgYXBwbHkoYm90bGVmdCwgMSwgZnVuY3Rpb24oeCl7cmV0dXJuKAogICAgb3dpbih4WzFdICsgYygwLCA1MCksIHhbMl0gKyBjKDAsIDUwKSkKICApfSkKKQpiZWlfaG9sZSA8LSBiZWlbY29tcGxlbWVudC5vd2luKGJlaV93aW4pXQpiZWlfc2FtcCA8LSBiZWlbYmVpX3dpbl0KYmVpX3dpbmRvd19mdWxsIDwtIFdpbmRvdyhiZWkpCgpwbG90KGJlaV9ob2xlLCBtYWluID0gJ09ic2VydmVkIFN1YnJlZ2lvbicsIHBjaCA9ICcuJywgY29scyA9ICdibGFjaycpCmBgYAoKYGBge3IgYmVpaG9sZW1lc2h9CmJlaV9ob2xlX21lc2ggPC0gaW5sYS5tZXNoLjJkKGNiaW5kKGJlaV9ob2xlJHgsIGJlaV9ob2xlJHkpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjdXRvZmYgPSA1MCwgbWF4LmVkZ2UgPSBjKDUwLCAxMDApLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBsb2MuZG9tYWluID0gYmVpX2RvbWFpbikKcGxvdChiZWlfaG9sZV9tZXNoLCBhc3AgPSAxKQpwbG90KFdpbmRvdyhiZWkpLCBib3JkZXIgPSAnYmx1ZScsIGFkZCA9IFRSVUUpCnBvaW50cyhiZWlfaG9sZSwgcGNoID0gJy4nLCBjb2wgPSAnYmx1ZScpCmBgYAoKYGBge3IgYmVpaG9sZWxnY3AsIGNhY2hlID0gVFJVRX0KYmVpX2hvbGVfc3BkZiA8LSBhcy5TcGF0aWFsUG9pbnRzLnBwcChiZWlfaG9sZSkKIyBDSEVDSyBQUklPUlMhCm1hdGVybl9ob2xlIDwtIGlubGEuc3BkZTIucGNtYXRlcm4oYmVpX2hvbGVfbWVzaCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwcmlvci5zaWdtYSA9IGMoMC4xLCAwLjk5KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwcmlvci5yYW5nZSA9IGMoNSwgMC4wMSkpCmNtcF9ob2xlIDwtIGNvb3JkaW5hdGVzIH4gbXlTbW9vdGgobWFwID0gY29vcmRpbmF0ZXMsIG1vZGVsID0gbWF0ZXJuX2hvbGUpICsgSW50ZXJjZXB0CmJlaV9ob2xlX2xnY3AgPC0gbGdjcChjbXBfaG9sZSwgYmVpX2hvbGVfc3BkZikKbGFtYmRhX2hvbGUgPC0gcHJlZGljdChiZWlfaG9sZV9sZ2NwLCBwaXhlbHMoYmVpX2hvbGVfbWVzaCksIH4gZXhwKG15U21vb3RoICsgSW50ZXJjZXB0KSkKCiMgUGxvdCBwb3N0ZXJpb3IgbWVhbnMgYW5kIHBvc3RlcmlvciBzZC4KcGxvdChsYW1iZGFfaG9sZSkKcGxvdChXaW5kb3coYmVpX2hvbGUpLCBib3JkZXIgPSAnd2hpdGUnLCBhZGQgPSBUUlVFKQpwb2ludHMoYmVpX2hvbGUsIHBjaCA9ICcuJywgY29sID0gJ3doaXRlJykKcGxvdChsYW1iZGFfaG9sZVsnc2QnXSkKcGxvdChXaW5kb3coYmVpX2hvbGUpLCBib3JkZXIgPSAnd2hpdGUnLCBhZGQgPSBUUlVFKQpwb2ludHMoYmVpX2hvbGUsIHBjaCA9ICcuJywgY29sID0gJ3doaXRlJykKYGBgCgpgYGB7ciBiZWlzYW1wfQpwbG90KGJlaV93aW5kb3dfZnVsbCwgbWFpbiA9ICdPYnNlcnZlZCBTYW1wbGUnKQpwbG90KGJlaV93aW4sIGFkZCA9IFRSVUUpCnBsb3QoYmVpX3NhbXAsIHBjaCA9ICcuJywgY29scyA9ICdibGFjaycsIGFkZCA9IFRSVUUpCmBgYAoKYGBge3IgYmVpc2FtcG1lc2h9CmJlaV9zYW1wX21lc2ggPC0gaW5sYS5tZXNoLjJkKGNiaW5kKGJlaV9zYW1wJHgsIGJlaV9zYW1wJHkpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjdXRvZmYgPSA1MCwgbWF4LmVkZ2UgPSBjKDUwLCAxMDApLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBsb2MuZG9tYWluID0gYmVpX2RvbWFpbikKcGxvdChiZWlfc2FtcF9tZXNoLCBhc3AgPSAxKQpwbG90KFdpbmRvdyhiZWkpLCBib3JkZXIgPSAnYmx1ZScsIGFkZCA9IFRSVUUpCnBsb3QoV2luZG93KGJlaV9zYW1wKSwgYm9yZGVyID0gJ2JsdWUnLCBhZGQgPSBUUlVFKQpwb2ludHMoYmVpX3NhbXAsIHBjaCA9ICcuJywgY29sID0gJ2JsdWUnKQpgYGAKCmBgYHtyIGJlaXNhbXBsZ2NwLCBjYWNoZSA9IFRSVUV9CmJlaV9zYW1wX3NwZGYgPC0gYXMuU3BhdGlhbFBvaW50cy5wcHAoYmVpX3NhbXApCiMgQ0hFQ0sgUFJJT1JTIQptYXRlcm5fc2FtcCA8LSBpbmxhLnNwZGUyLnBjbWF0ZXJuKGJlaV9zYW1wX21lc2gsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcHJpb3Iuc2lnbWEgPSBjKDAuMSwgMC45OSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcHJpb3IucmFuZ2UgPSBjKDUsIDAuMDEpKQpjbXBfc2FtcCA8LSBjb29yZGluYXRlcyB+IG15U21vb3RoKG1hcCA9IGNvb3JkaW5hdGVzLCBtb2RlbCA9IG1hdGVybl9zYW1wKSArIEludGVyY2VwdApiZWlfc2FtcF9sZ2NwIDwtIGxnY3AoY21wX3NhbXAsIGJlaV9zYW1wX3NwZGYpCmxhbWJkYV9zYW1wIDwtIHByZWRpY3QoYmVpX3NhbXBfbGdjcCwgcGl4ZWxzKGJlaV9zYW1wX21lc2gpLCB+IGV4cChteVNtb290aCArIEludGVyY2VwdCkpCgojIFBsb3QgcG9zdGVyaW9yIG1lYW5zIGFuZCBwb3N0ZXJpb3Igc2QuCnBsb3QobGFtYmRhX3NhbXApCnBsb3QoV2luZG93KGJlaSksIGJvcmRlciA9ICd3aGl0ZScsIGFkZCA9IFRSVUUpCnBsb3QoV2luZG93KGJlaV9zYW1wKSwgYm9yZGVyID0gJ3doaXRlJywgYWRkID0gVFJVRSkKcG9pbnRzKGJlaV9zYW1wLCBwY2ggPSAnLicsIGNvbCA9ICd3aGl0ZScpCnBsb3QobGFtYmRhX3NhbXBbJ3NkJ10pCnBsb3QoV2luZG93KGJlaSksIGJvcmRlciA9ICd3aGl0ZScsIGFkZCA9IFRSVUUpCnBsb3QoV2luZG93KGJlaV9zYW1wKSwgYm9yZGVyID0gJ3doaXRlJywgYWRkID0gVFJVRSkKcG9pbnRzKGJlaV9zYW1wLCBwY2ggPSAnLicsIGNvbCA9ICd3aGl0ZScpCmBgYAoKCiMgUmVmZXJlbmNlcwoK